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Abstract 


Whereas diffusion NMR can probe the structural configurations of microscopic environments in biological 
tissue, relaxation can provide complementary information on their chemical composition. This chapter 
considers experiments in which diffusion and relaxation properties are sampled simultaneously by varying 
multiple acquisition parameters. As such, correlations between the diffusion and relaxation can be established, 
providing an altogether more complete picture of heterogeneous tissue. 


7.1. Introduction 


The previous chapters have dealt with the use of magnetic resonance to track the translational motion of 
molecules. Multi-dimensional diffusion MR encodes the signal in multiple dimensions; e.g., two successive 
pulsed-gradient pairs in different directions. If molecules experience hindered diffusion in locally anisotropic 
micro-domains, the signal attenuation will depend on the relative orientation of the pulses. As shown in Chapter 
3, this can thus be used to characterise local anisotropy even if the orientation distribution of the domains is 
globally random. 


The experiment with two successive pulsed-gradient pairs in different directions is part of a much broader class 
of experiments based on the concept of multiplexing! in which multiple signals are combined into one signal. In 
the context of MR, this means that the signal is encoded as a function of multiple experimental parameters. In 
contrast to separately acquiring each dimension, multi-dimensional MR can be used to estimate the joint 
distribution and as such establish correlations between the dimensions. 
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2 Advanced Diffusion Encoding Methods in MRI 


Peemoeller et al. (1981)2 describes one of the earliest multi-dimensional MR experiments of relaxation times in 
wet hen egg-white lysozyme. They analysed the data using a graphical method called curve peeling, which does 
not require non-linear curve fitting to extract exponentials. In multi-dimensional experiments, a 
multiexponential decomposition can then be graphically performed in each dimension. English et al. (1991) 
later introduced numerical methods, which were not biased by a priori assumptions on the number of 
components and yielded reproducible results. Since then, many approaches have been proposed to numerically 
estimate correlation spectra”, 


This chapter focuses on the correlation of diffusion and relaxation, for which the mathematical formulation was 
already presented in the 70s8. Heink et al. (1991)? observed an interaction between diffusion and relaxation: the 
estimated apparent diffusivity turned out to depend on the echo time (TE), and this was initially treated as an 
artifact. Van Dusschoten et al. (1995) 10 later exploited this as a contrast to sensitise the signal to the different 
physical phenomena underlying diffusion and relaxation. Since then, diffusion-relaxation correlation has been 
applied to a wide range of porous media, e.g. in sedimentary rocks11-14, cheeses and other dairy products!>!6, 
liquid crystals” !”, glass beads18,19 drugs”, crude oils*!22, hydrogels23, plants!®2425, and yeast cells26_ It has 
even been implemented to be used with inhomogeneous low magnetic fields without pulsed gradients: a 
permanent magnet with a static gradient and a basic NMR spectrometer greatly facilitates application in large 
substances!®, 


Multi-dimensional diffusion-relaxation correlation experiments are becoming more popular in human 
MRI27-38_ In current clinical practise, routinely acquired MR-images still largely show qualitative tissue- 
contrast, with intensities arbitrarily scaled according to different MR-phenomena. In contrast, quantitative MRI 
aims to extract reproducible measures more directly related to tissue properties. Diffusion-relaxation correlation 
increases the quantitative potential of MRI and as such can have important clinical ramifications in disease 
management and diagnosis. This chapter will focus on diffusion-relaxation correlation in animal and human 
tissue, in the light of its potential implications on clinical MRI. 


The chapter is structured as follows. Section 7.2 will briefly answer the question “what is relaxation?”, 
complementing the extensive discussions in the previous Chapters on the process of diffusion. Section 7.3 will 
give an overview of studies that have performed diffusion-relaxation correlation in animal tissues, answering the 
question “why complement diffusion with relaxation?” This Section is structured according to four commonly 
reported key benefits of diffusion-relaxation correlation over either diffusion- or relaxation-MR alone: 1) it 
provides additional voxel-wise information, 2) multiple components per voxel can be better disentangled, 3) it 
can stabilise the fitting in either the diffusion or relaxation dimension, and 4) it can ameliorate the issue of 
diffusion MR providing T2-weighted signal fractions in the case of multiple components. Section 7.4 briefly 
discusses “how to measure diffusion and relaxation simultaneously?”, where the reader is also referred to 
Chapter 9 for a more detailed description on acquisition. Section 7.5 considers the topic “how to estimate 
diffusion and relaxation?” introducing free inversion and model-based estimation (Chapters 11 and 12). Finally, 
Section 7.6 will review the question “how to validate diffusion-relaxation correlation?” and Section 7.7 discusses 
limitations and future outlooks. 


7.2. What is relaxation? 


Relaxation is the process of spins returning to their thermal equilibrium state with time, after they have been 
brought into an excited state. T1 and T2 relaxation describe two different processes affecting the nuclear 
magnetisation; both have a complex relation to material properties. Relaxometry has been extensively used to 
quantitatively characterise a wide range of materials (see e.g.3? for an overview in liquids and*° in brain). In 
biological tissues, relaxation can originate from phenomena on different scales: the molecular nanometer scale, 
the cellular micrometer scale, and the MRI resolution millimetre scale#!~49. Here, we give a brief overview of 
the different relaxation processes. 
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7.2.1. Longitudinal relaxation 


Longitudinal relaxation, or spin-lattice/thermal/T1 relaxation, is the recovery of the longitudinal (i.e. in z- 
direction) nuclear spin magnetisation towards the thermal equilibrium. The T1 value is the time it costs for the 
longitudinal magnetisation to recover to 63% (i.e. 1 — exp(—1)) of its initial value after a 90° RF pulse. On the 
molecular scale, longitudinal relaxation involves a redistribution of spin-up and spin-down populations for 
which the energy gained from the RF pulse must leave the spin system. To this end, the spins have to exchange 
energy with their surroundings (the lattice) as spontaneous emission is unlikely in the frequency range of NMR. 
The most efficient longitudinal relaxation occurs when protons or electrons on the same or neighbouring 
molecules cause locally fluctuating magnetic fields near the Larmor frequency“. The optimal ‘interaction 
distance’ in MRI is in the order of nanometers, and the values of T1 are mainly given by cross-relaxation 


between water and the macromolecules in tissue*. 


7.2.2. Transverse relaxation 


Transverse relaxation, or spin-spin relaxation, is the decay of the net transverse magnetisation (i.e. in the xy- 
plane) due to dephasing of the individual spins. The T2 value is the time it costs for the transverse relaxation to 
decay to 37% (i.e. exp(—1)) of its initial value after a 90° RF pulse. The T2 effect, in contrast to T1, can originate 
from different scales and is mostly caused by a variation in the magnetic field. On the molecular scale, the 
magnetic field changes due to random rotation of the molecules (and the associated dipole fields of their 
protons). In the presence of large molecules, this motion slows down and the relaxation rate (Ro = 1/ T>) 
increases**4. In addition, if a spin that contributes to the net transverse magnetisation exchanges energy with the 
lattice following T1 relaxation, then it will also lose phase coherence with other spins. 


Magnetic susceptibility is an inherent property of a material and describes the relation between an externally 
applied magnetic field and the field within the object. In composite materials the magnetic susceptibility varies 
with position in space, which leads to spatial variation of the local magnetic field. On the cellular scale (also 
referred to as microstructural or mesoscopic scale), the presence of paramagnetic substances induces an 
additional magnetic field that varies in space and time. This causes additional dephasing of spins not only in the 
immediate neighbourhood of the substance, but also on larger length scales depending on the size and geometry 
of the compartment in which the substance is residing. In porous media such as sedimentary rocks, 
paramagnetic ions on or near the surface cause enhanced relaxation of fluid molecules near the surface, also 
called surface relaxation*®. This has been used to study the surface-to-volume ratios of pores*’. In biological 
tissue, iron associated with hemoglobin in blood or in cells (e.g. in deep gray matter nuclei in the brain) cause 
additional transverse relaxation. 


Even if the field inhomogeneities are static (that is, they don’t change over time), the spins probe different 
magnetic susceptibilities over time if they exhibit non-negligible diffusion. As such, the measurement becomes 
sensitive to the spatial correlations of the susceptibility. Part of this relaxation cannot be fully refocused, as the 
effect of diffusion cannot be undone. The behaviour of the transverse relaxation as a result of diffusion in 
inhomogeneous fields has been investigated in several works!®-48->0, 


Finally, on the macroscopic scale of the MRI voxel size, inhomogeneities can be caused by differences in 
macroscopic susceptibility between tissues, e.g. at the air-tissue interface near sinuses in the brain, by poor 
shimming, and by the sample shape. The effect of diffusion can be neglected at this scale, and dephasing due to 
these inhomogeneities can in principle be fully rephased by the spin-echo experiment. The classic equation for 
NMR in chemical solutions separates relaxation from macroscopic origins (denoted by R3): 


R¥ = R +R, [2] 
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7.3. Why complement diffusion with relaxation? 


7.3.1. Additional information 


From the previous section, it is obvious that relaxation provides a rich source of information that is 
complementary to diffusion; whereas diffusion NMR provides information on molecular motion, relaxation 
provides information on the molecular interactions of the system of spins, amongst others. Hutter et al. (2018 
showed scatter plots of estimated voxel-wise T1, T2*, and diffusivity in the living human brain, revealing the 
independence of these tissue parameters and thus the benefits of acquiring multiple dimensions (Figure 7.1a). 
Simultaneous quantification of voxel-wise NMR properties (the NMR fingerprint?!) can improve the 
classification of tissues, and as such the diagnostic capability in medical applications. This was for example 
demonstrated by Yu et al. (2017)°2, who performed T1, T2, and diffusivity mapping in prostate and showed that 
tissue manually classified as prostate cancer had lower estimated T1, T2, and D than normal tissue. Combined 
contrasts (specifically the estimation of voxel-wise D and T2, Figure 7.1b) had a better separation between 
different tumour grades than single contrasts. 


)28 


Figure 7.1. Diffusion-relaxation correlation gives additional information. Here, one D and one T2 and/or T1 
per voxel were estimated and then distributions across voxels are shown. 
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a). In healthy controls, scatter plots between voxel-wise T2, T1, and D reveal independence between parameters. 
Copyright 2018 


Adapted from?® under the Creative Commons license (http://creativecommons.org/licenses/by/4.0/) 
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b). In prostate cancer patients, D and T2 were significant predictors for differentiating between histologically proved high- or 
intermediate-grade tumours and low-grade tumours. 


Adapted from”? with permission from RSNA, Copyright 2017. 


7.3.2. Resolving multiple components 


For simple systems such as water in aqueous solutions, all the molecules sense the same environment and the 
material can thus be described with a single T1, T2, and D (also the assumption in most MR fingerprinting 
approaches). For complex systems such as porous rock, food, and biological tissue, the molecules can reside in 
different microscopic environments (e.g. inside and outside cells in tissue). As such, different molecules can 
sense different microscopic environments, and the material is better described by distributions of T1, T2, and D 
within a voxel (as opposed to Figure 7.1 where distributions of voxel-wise T1, T2, and D across voxels are 
shown). 


The 1D distributions of NMR observables in heterogeneous materials are often quite broad, and hence 1D 
experiments might struggle to resolve separate components in these materials. Multi-dimensional experiments, 
in contrast, can rely on measurements of different underlying physical processes (e.g. relaxation and self- 
diffusion) to disentangle compartments and thus have the potential to improve the resolving power of NMR. 
Instead of just acquiring different 1D experiments, multi-dimensional correlation experiments - in which the 
signal attenuation is acquired by varying multiple experimental parameters simultaneously - are required to 
resolve the joint distribution of the different NMR observables. Indeed, this principle was first demonstrated in a 
2D T1-T2 experiment in excised rat muscle*, where independent 1D experiments failed to resolve all 
components. 


Diffusion-relaxation correlation, which is the focus of the current Chapter, has shown to be a powerful tool to 
disentangle components that are different in size/shape and chemical composition. The next subsections will 
focus on experiments to investigate compartmentalisation in animal and human tissue. 
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7.3.2.1. In vitro NMR 


Nervous tissue — which is commonly compartmentalised into water residing between the myelin sheets, within 
the axons, and outside the axons (Figure 7.2a) — has been studied extensively with diffusion-relaxation 
correlation. Beaulieu et al. (1998)°> sought to provide evidence of compartmentalisation in nerve water T2 and 
D by using three different excised nerves of the garfish; non-myelinated olfactory nerve, trigeminal nerve 
myelinated by Schwann cells, and optic nerve myelinated by oligodendrocytes. They indeed found three T2 
components in myelinated nerves but only two in the non-myelinated nerves. They did not find any dependence 
of the measured apparent diffusivity and anisotropy on the T2 component of origin. Stanisz and Henkelman 
(1998)>4, in contrast, found that diffusion was more anisotropic in the longest T2-component than in the 
shortest T2-component in the bovine optic nerve, where the latter is often associated with water confined within 
myelin. Peled et al. (1999)°° found three T2 compartments in the frog peripheral nerve, where the shortest T2 
compartment showed some restriction, the intermediate T2 compartment showed relatively high and 
unrestricted diffusion, and the longest T2 compartment showed restricted diffusion. The latter two 
compartments were assigned to extracellular and intracellular (most likely axonal) water. Andrews et al. 
(2006)°° studied the diffusion-T2 characteristics of myelin water specifically, where they also used T1 
characteristics to nullify non-myelin water. They found a lower perpendicular diffusivity, and thus higher degree 
of anisotropy, than Stanisz and Henkelman (1998)>4, 


Other excised tissues and cells have been studied with diffusion-relaxation correlation NMR. Pfeuffer et al. 
(1998)°7 used diffusion-T2 relaxometry to investigate the intra and extracellular water signal and water 
permeability in Xenopus Laevis oocytes and glial cells. Li et al. (1998)°° report that integrated D-T2 experiments 
result in better parameter determinations in a three-pool model of water in human blood. Seland et al. (2005)? 
show two components in rat myocardium with distinctly different diffusivity and T1, which they assign to intra- 
and extra-cellular water. Mailhiot et al. (2018)®9 show two components in porcine articular cartilage with a 
common D value but different T2 values. 


7.3.2.2. Animal MRI 


When combined with imaging, one can spatially map different components to study variations across the 
tissue!®. Several studies have performed diffusion-relaxation correlation MRI in animals, both ex vivo and in 
vivo. 


Ex vivo diffusion-relaxation correlation MRI experiments have been reported in mouse and rat spinal cord, and 
rat brain, amongst others. Kim et al. (2017)°! measured diffusion and T2 relaxation in ex vivo mouse spinal 
cords, which revealed two different peaks to disentangle white matter, which contains the axons, and gray 
matter, which contains mostly cell bodies. They found an additional peak in injured cords that could not be 
found in healthy controls. Compared to the 2D approach, there was considerable ambiguity in the 1D diffusion 
and 1D relaxation spectra which made it difficult to interpret them in a meaningful way. Benjamini and Basser 
(2017) performed D-T2 experiments in dissected rat spinal cord, and found evidence for three distinct 
components within the gray matter, which they assign to intracellular/neuronal soma (lower D), interstitial 
(higher D), and myelin-associated (lower T2) spaces. They show that the intracellular T2 is slightly longer than 
the extracellular T2. In white matter, two components were found that were assigned to interstitial and myelin- 
associated spaces, and an additional two components that were assigned to intracellular spaces: one longer T2- 
component that was also found in gray matter and assigned to glial soma water, and one shorter T2-component 
assigned to intra-axonal water. De Santis et al. (2016)® correlated diffusion tensors and T1-values in two excised 
porcine nerve fibres with very different myelination that are crossing each other at an angle. Both fibres exhibit 
very different T1 values in the crossing area when combined diffusion-relaxation information is into account, 
whereas 1D T1-relaxometry fails to recover two distinct T1s. This was also seen in perfused rat brain, where the 
estimated T1s were different in two crossing fibres (cingulum and genu of the corpus callosum). 
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In vivo diffusion-relaxation correlation MRI experiments have been reported in rat brain, rat trigeminal nerve, 
and cat brain, amongst others. De Santis et al. (2016)? extended their D-T1 MRI experiment to in vivo rat 
brains, and showed that the cingulum consistently exhibited higher T1 than the corpus callosum and that the 
more restricted population had a higher T1 in the cortex and the gray matter, but a smaller T1 in the white 
matter. van Dusschoten et al. (1996)54 performed diffusion-T2 correlation in healthy cats and cats with ischemic 
regions following a middle cerebral artery occlusion. They found that the properties of cerebral spinal tissue and 
brain tissue could be more reliably analysed than with separate measurements. In addition, whereas they did 
find a lower diffusivity for the compartment associated to the ischemic region, they did not find T2 changes in 
the compartment associated to ‘normal brain tissue in ischemic regions. Does and Gore (2000)® imaged 
healthy and ischemic rat trigeminal nerve and brain, and found no proof of compartmentalisation in brain, 
possibly due to relatively rapid exchange between intra- and extra-cellular compartments. Two compartments 
were found in trigeminal nerve: water with longer T2 was found to be more restricted in perpendicular direction 
than water with shorter T2, possibly corresponding to extracellular and intracellular water. Myelin water was off- 
limits with the acquisition protocol used. Qin et al. (2009)®° show that simulations with a longer T2 for the 
extracellular compartment support the observed behaviour of TE dependence in rhesus monkey. In contrast, 
results of De Santis et al. (2016)®” in rat brain support a longer T2 for intra-axonal water. 


7.3.2.3. Human MRI 


Diffusion-relaxation correlation has also been performed in the living human brain and other organs. One major 
complication of in vivo human imaging is the long scan time; it takes significantly more time to adequately 
sample a 2D or 3D acquisition space compared to 1D (see also Section 7.4.3). Mulkern et al. (2000)68 performed 
T1-D experiments to make separate measurements of T1 relaxation of the components with slow- and fast 
diffusion in the human brain (often associated to intra- and extra-cellular components). They did not find a 
significant difference, which led them to conclude that either the T1’s within the compartments are similar, or 
the rate of water exchange is fast enough to mix the T1-relaxation behaviour of the compartments but slow 
enough to show a diffusion-behaviour consistent with multiple compartments. De Santis et al. (2016)32 applied 
their diffusion-tensor T1 experiment to human brain, reproducing the result of different T1 values in crossing 
tracts. Moreover, they found a lower inter-subject variability for parameters estimated with T1-D correlation 
compared to T1-relaxometry alone. 


Using D-T2 correlation, de Almeida Martins et al. (2019)® resolve five-dimensional diffusion-relaxation 
distributions in the living human brain (a distribution of T2s and axially symmetric tensors per voxel, the latter 
including perpendicular and parallel diffusivities, and two angles defining the first eigenvector). This high 
dimensionality allows them to resolve multiple microscopic environments within each voxel, e.g. in voxels with 
partial volumeing of CSF, gray matter, and white matter (Figure 7.2c). Regarding the separation of intra- and 
extra-axonal T2 using D-T2 correlation, Lin et al. (2018)”9 report that simulations of two hypotheses support 
their observations in the living human brain equally well: either the intra-axonal T2 is larger than the extra- 
axonal T2, or both spaces have intrinsically similar T2 but exchange with the myelin compartment causes a 
measurable difference. They hypothesise that the axonal membrane reduces the exchange rate between myelin 
water and intra-axonal water, whereas the larger interface with the extra-axonal increases the exchange with the 
extra-axonal water, causing the extra-axonal T2 to be reduced to a larger degree. Peled et al. (1999)>» attribute 
the lower extracellular T2 to simple chemical composition effects; a higher collagen concentration in the 
extracellular space drives water relaxation via chemical exchange. Several studies27;31,35 which correlate T2- 
relaxation and diffusion tensors, support the hypothesis of a larger intra-axonal T2 (Figure 7.2d). Various 
works>!)71,72 furthermore report a dependence of the voxel-wise- or intra-axonal T2 on the angle of the 
underlying tract relative to the main magnetic field. Ning et al. (2019)73 show that in crossing fibre areas, 
different bundles were estimated to have different relaxation rates, which may have important implications for 
diffusion tractography. Lemberskiy et al. (2018)’4 perform D-T2 experiments in prostate and use the difference 
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in T2 of the cellular and luminal compartments to disentangle their contributions to the time-dependent 
diffusion coefficient. 


D-T2* experiments have also been performed in humans. Slator et al. (2019)75 studied the human placenta and 
found three components with clearly different diffusivities but similar T2*. These compartmental properties 
changed in patients with placentally-mediated pregnancy complications. Kleban et al. (2019)7° aim to 
disentangle intra- and extra-axonal T2* using diffusion weighting, and find a significantly longer T2* for the 
intra-axonal compartment when the underlying fibre is perpendicular to the main magnetic field. In addition, 
they characterise compartmental variations in magnetic susceptibility by estimating the susceptibility-induced 
frequency offset, which can be interpreted as a proxy for the anisotropic susceptibility of the myelin sheath. In 
1D T2*-relaxometry, including the frequency offset has shown to improve the estimation of the myelin water 
fraction’’, but it has also been argued that this comes at the cost of increased complexity of the fitting procedure 
and may therefore be ignored at field strengths of 3T and lower’®, Joint diffusion-T2* relaxometry can provide 
new avenue? to characterise the complementary microscopic information in the estimated frequency offset more 
reliably” °. 


Figure 7.2. Diffusion-relaxometry can better disentangle different compartments. An example in human 
brain white matter and cerebrospinal fluid (CSF). 


a). Schematic compartmentalisation of nervous tissue into myelin water, extra-axonal water, and intra-axonal water. Myelin is a 
stack of lipid bilayers and forms a concentric wrapping around the axon to ‘insulate’ the axon for the transmission of electrical signals; 
water is contained in the regions of the hydrophilic headgroups. 
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b). Schematic D-T2 spectrum for white matter. 
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c). Left: diffusion-relaxation distributions shown as scatter plots in a 3D space of logarithms of the transverse relaxation rate R2, 
isotropic diffusivity Diso, and axial-radial diffusivity ratio D|// D | The diffusion tensor orientation (q,f) is color-coded as [R,G,B] 
= [cos ó sin 9, sin $ cos 0, cos 0]|D|| — D ( /max(D||,D 1) and the circle size is proportional to the weight of the component. The 
distributions are subsequently divided into ‘bins’ named ‘Big; “Thin, and “Thick that loosely capture the diffusion features of 
cerebrospinal fluid CSF, white matter WM, and gray matter GM, respectively. Right: Parameter maps of average per-bin means (color) 
of transverse relaxation rate, isotropic diffusivity, squared anisotropy and diffusion tensor orientation, color-coded as [R,G,B] = [Dxx, 
Dyy, Dzz]/max(Dxx, Dyy, Dzz), where Dii are the diagonal elements of laboratory-framed average diffusion tensors estimated from the 
various distribution bins. Brightness indicates the signal fractions corresponding to the ‘Big’ (row 1), “Thin (row 2), and “Thick (row 3) 
bins. 

Copyright 2019 

Adapted from®? under the Creative Commons license (http://creativecommons.org/licenses/by/4.0/) 
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d). There is some debate in the literature, but recent results show that T2 is longer in the more-restricted intra-axonal space in the 
living human brain. 


Adapted from27 with permission from Elsevier, Copyright 2018. 
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7.3.3. Improve single-contrast parameter estimation 


Reliable estimation of parameters is challenging with the typically low SNR in diffusion NMR and MRI. In order 
to obtain a good contrast, sufficiently strong diffusion weighting and/or long diffusion time has to be applied, 
resulting in a decay of the signal and thus SNR. Adding measurements along an independent dimension not only 
increases the effective SNR®*79, but can also improves the diffusion estimates. 


van Dusschoten et al. (1995)!9 used the difference in T2 relaxation times in apple parenchyma and mung bean 
seedlings to discriminate diffusion constants. They observed that the errors of the estimated diffusion constants 
were much smaller than when a multi-component model was fitted to the diffusion data alone. Celik et al. 
(2013)®° show that the accuracy of estimates and their stability with respect to noise is markedly improved in 
general when a second independent dimension is introduced (Figure 7.3a). While this is showcased for T1-T2 
correlation, the results are generalisable to diffusion-relaxation correlation. 


In the living human brain it has turned out to be challenging to achieve a unique, biophysically plausible 
solution for a simple two-component model (intra- and extra-axonal space) using diffusion MRI data alone. 
Specifically, an analysis of this model revealed a degeneracy in parameter estimation on data up to moderate 
diffusion weightings; two biophysically plausible solutions were shown to be equally supported by the data81, 
Several studies?”:>> show that by measuring both diffusion and T2-relaxation, the two compartments can be 
disentangled with increased confidence, where high b-values further improve the result (Figure 7.3b). 


Figure 7.3. Diffusion-relaxometry can improve single-contrast parameter fitting. 


weights 
° _ ° 
—— 
u 


T, (ms) 


7 
T2 


T, (ms 


T, Projection (ms) 


a). Results in a two-component gel, using nonparametric inversion of T1-T2 data (see Section 7.5.1.1). (i) T2 distributions for 
individual gels combined into a single histogram. (ii) T2 histogram for the two-component sample via inversion, revealing only a single 
compartment. (iii) T1-T2 histogram for the two-component sample via inversion. (d) T2 projection of the T1-T2 histogram, now 
revealing both compartments. 


Figure adapted from®? with permission from Elsevier, Copyright 2013. 
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b). Results in a two-component D-T2 simulation of white matter, using model-based estimation (Section 7.5.1.2). The ground 
truth is indicated by the black dot, and the left and right panel show T2 vs the perpendicular and parallel diffusivity of the modelled 
tensor, respectively. Including higher b-values reduces the observation of a ‘second solution branch (bottom). T2 in ms, D in mm2/s. 


7.3.4. Weighted signal fractions 


Using diffusion NMR alone, multiple studies report signal decays in brain consistent with the presence of intra- 
vs extra-cellular compartments. However, most of these studies report a discrepancy between the estimated 
signal fractions and literature values for the volume fractions corresponding to these compartments. Niendorf et 
al. (1996)®2 argue that a difference in T2 between the intra- and extra-cellular space makes the estimated signal 
fractions weighted according to their T2, which could cause this discrepancy. Lampinen et al. (2019)34 report 
this as a major challenge to estimate the neurite volume fraction, i.e. the fractional voxel volume occupied by 
neurites. Several studies that combine diffusion and T2 measurements find values that are in better agreement 
with literature values, e.g.27©2-83, However, even if the signal fractions were corrected for relaxation, converting 
them to volume fractions requires knowledge of the number of detected spins per unit volume in each of the 
compartments, and this property can be very different for, e.g., myelin, intra- and extra-axonal spaces, and CSF. 
Claims of measured neurite densities or volume fractions with MRI should therefore be interpreted with caution. 


7.4. How to measure diffusion and relaxation simultaneously? 


MRI experiments generally consist of an excitation block (commonly a 90° radiofrequency (RF) pulse), a 
weighting block, and an acquisition/readout block. Here, the aim is to give some intuition as to how different 
weightings can be used to create diffusion-relaxation correlation experiments, and we refer to Chapter 9 for a 
more extensive discussion on weightings and readouts. In diffusion MRI, the spin echo (SE) sequence is the most 
commonly used; diffusion gradients are placed on both sides of the 180° pulse to sensitise the acquisition to 
diffusion (Fig 484). The SE sequence is simple and elegant, and we will use this as an example of how diffusion 
and relaxation can be measured simultaneously. 
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7.4.1. Measuring diffusion and relaxation with the spin-echo sequence 


7.4.1.1. A brief recap of the diffusion-gradient spin-echo sequence 


The 90° RF pulse will bring the net magnetisation vector in the transverse plane in a given slice, with a fraction 
of the spins being in phase. From there, the transverse magnetisation component starts precessing around the z- 
axis with the Larmor frequency. We can write the transverse magnetisation as a complex number that varies over 
time, with the x-component (sinusoidal) the real and the y-component (cosinusoidal) the imaginary part. The 
longitudinal magnetisation is zero just after the 90° RF pulse. Due to transverse relaxation effects at different 
scales (Section 7.2.2), each individual spin will precess with a slightly different frequency, and the dephasing 
causes the magnitude of the transverse magnetisation to decay. 


Then, diffusion weighting is applied by switching on gradients in different directions. Here, a pulsed-gradient is 
illustrated (the shaded rectangle in Figure 7.4). The gradient causes the total magnetic field to vary spatially: 
spins precess a bit slower or faster depending on their position. At the end of the gradient waveform, the spins 
have a position-dependent phase shift (they ‘dephased; where we ignore relaxation). After the application of the 
diffusion gradient, all spins start to precess again with their ‘normal’ frequency. 


At time point TE/2 (with TE the echo time) a 180° pulse is applied, ‘flipping’ the phase of the spins. This pulse is 
part of the spin-echo sequence and has now started the correction for reversible transverse dephasing due to 
chemical shift offset and BO inhomogeneity. This pulse will yield an echo at time TE, with the signal decay now 
solely originating from irreversible transverse relaxation. Note that it also inverts the sign of the phase shift 
caused by the diffusion gradient. The fraction of the longitudinal magnetisation that has grown back will be 
flipped along the negative axis. 


A second diffusion gradient waveform is applied which now tries to undo the phase shift caused by the first 
gradient. This is the essence of the diffusion encoding: when the molecules remain at the same location the 
effects of both gradients indeed cancel out, but when these molecules have diffused the effect of the first gradient 
cannot be undone. The net dephasing of a spin is dependent on its displacement in the direction of the diffusion 
gradient, and the strength and duration of the gradients. This dephasing of spins will result in a signal drop in 
addition to the decay from transverse relaxation at time point TE. Just after the echo, the signal will decay again. 
The longitudinal magnetisation will continue growing until the next 90° RF pulse at the repetition time (TR). 


The following sections discuss how relaxation can simultaneously be measured by varying multiple parameters. 
Here, the classical description as in Eq. 2 will be adopted. 


7.4.1.2. Transverse relaxation 


From the diffusion-weighted SE sequence it quickly becomes evident that varying the TE will give different T2 
weightings to the signal. In the case of a single compartment and fully recovered magnetisation just before the 
90° RF pulse (e.g. during the first TR, denoted by So), the signal at time TE resulting from transverse relaxation 
(ignoring diffusion) is 


S(TE) = Svexp(-75)- [3] 


In the case of different compartments, the compartment with the longest T2 will have the largest signal fraction. 
Various works27;35,57,61 use pulsed gradient (PG, or Stejskal-Tanner encoding84) SE, and vary the TE. When 
additionally varying the gradient strength and direction, the resulting data will be 4-dimensional in acquisition 
space (TE,b,0,¢). Others?3-85 additionally vary the waveforms, to generate tensor-valued diffusion encoding 
(Chapter 3). They limit the acquisitions to axially symmetric b-tensors, with b4 = (b|| _ b | )/ (b|| —2b,) 
characterising the anisotropy. As such, the acquisition space becomes 5-dimensional (TE,b,b,,6,¢). 
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Figure 7.4. The spin-echo sequence can be used for diffusion-relaxation correlation. GR, Gpg and Gs are the readout-, phase- 
encoding-, and slice-select gradients, respectively. ó and A characterise the pulse duration and separation, and G the gradient strength. 
The coloured arrows schematically represent different the different fields that spins experience due to the gradient being non-zero. 


The reversible spin dephasing is completely refocussed in a SE sequence at time TE. Boxerman et al. (1995) and 
Stables et al. (1998)8®87 proposed to shift the centre of the acquisition window in time relative to the centre of 
the spin echo by delay time TD, to achieve an ‘asymmetric’ SE sequence (ASE). In this way, the reversible 
dephasing contributes to the signal decay. Kleban et al. (2019)’© combine diffusion PG with ASE to disentangle 
compartments according to their T2* and frequency difference. 


7.4.2.3. Longitudinal relaxation 


The T1 weighting can be modulated by varying TR and TE. In the steady state, which occurs during the second 
TR, the signal at time TE resulting from longitudinal relaxation (again ignoring diffusion) is 


S(TE,TR) = sdi = exp- EE) + exp| Tr] r) [4] 


If TE < TR, then TR — TE = TR, and Eq. 4 can be further simplified showing that the T1 modulation depends 
primarily on TR. This strategy to measure T1 is also referred to as saturation recovery (SR) SE; Seland et al. 
(2005)°? perform D-T1 correlation experiments using SR-PGSE. 


7.4.2. Measuring diffusion-relaxation with variations on the spin-echo 
sequence and other sequences 


Instead of a single 180° pulse and readout per excitation, the Carr-Purcell-Meiboom-Gill (CPMG) sequence 
applies multiple 180° pulses resulting in multiple echoes. CPMG can greatly speed up acquisition time and can 
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be placed before the diffusion weighting>’, e.g. to select the T2 range of the water proton magnetisation entering 
the diffusion measuring sequence, or after diffusion weighting®*. However, care should be taken regarding 
spurious echoes and imperfect refocusing. van Dusschoten et al. (1995)°8 use crushers to prevent spurious 
echoes, whereas Stanisz and Henkelman (1998)°4 use phase cycling. Several works>+88,89 use short 180 pulse 
intervals to minimise the effects of diffusion in background gradients (see also Section 7.6.5). 


Rather than varying TR to achieve different T1 weightings, Mulkern et al. (2000)®8 incorporate an inversion 
pulse prior to the 90° pulse in the diffusion-weighted SE sequence for simultaneous D-T1 measurement. The 
signal attenuation at time TE resulting from relaxation can in that case be written as 
S(TE,TR,TI) = Sl — 2exp(- rl + exp(- 71) ex0(- >) [5] 
a T1 T1 T2 J 
with TI the inversion time. 


Instead of a single 180° pulse in the SE experiment, a stimulated echo (STE) experiment uses two 90° pulses with 
a mixing time TM between the pulses; the advantage being that during TM the signal decays with (much slower) 
T1%°. This provides avenues to increase the time for diffusion weighting without losing too much signal because 
of T2, and to inherently add T1 weighting to the signal. However, crushers are needed to eliminate spurious 
echoes resulting in an SNR penalty. In addition, complications with cross relaxation (magnetization transfer) 
give an artefactual time-dependence of D for STE sequences when applied to water in biomaterials like brain 
tissue#!-45.91_ Diffusion-relaxation correlation using STE is adopted in several works! }89, 


7.4.3. Sampling the diffusion-relaxation acquisition space 


The strategy of sampling the MRI acquisition parameters can have a big impact on the ability to resolve the 
relevant information from the data. Relaxation and diffusion phenomena are frequently described by 
exponential decays (e.g. Eqs. 3 and 4), and the influence of the sampling extent and rate on exponential analyses 
has been studied for decades”. To optimally disentangle different decay rates (e.g. associated to different 
compartments), the signal should be recorded until it decays into the noise. While a linear spacing of data points 
seems the most straightforward, other - perhaps less intuitive — sampling schemes have been shown to be able to 
improve results. Bertero et al. (1984)?> show that with only 5 ‘geometrically spaced’ data points (i.e. t; = tid} 1i 
= 1,2,...,N), the restoration of up to 4 exponential components is improved compared to 32 linearly spaced data 
points (i.e. t; = tıd(i — 1),i = 1,2,...,N). 


For relaxometry, Weiss et al. (1980)? have studied optimal sampling schemes for T1 measurement, and Jones et 
al. (1996)?> have applied Cramer-Rao theory to find an optimal set of TEs that minimise the variance of the 
estimated T2 in case of a single compartment. They find that placing 22% of the sampling points at zero TE and 
78% at 1.28 T2 yields the best results. Beaulieu et al. (1998)23 resample linearly spaced echoes to approximate 
logarithmic sampling. Even though such non-linear echo spacing has been shown to improve accuracy of T2 
estimation, Does and Gore (2000)° point out that this may not be appropriate in many situations; imperfect 
refocusing pulses cause progressive loss of transverse magnetisation and can alternate relative T2 fractions or 
introduce non-monoexponential T2 behaviour. For diffusion, methods for optimal angular sampling have been 
introduced”, as well as strategies for optimal angular and radial sampling (i.e., in the case of different b- 
values)?®. More recently, such optimisations have also been investigated for b-tensor encoding?” 100, 


However, in the case of diffusion-relaxation correlation, the dimensionality of the MRI-acquisition space further 
increases and the number of combinations of MRI parameters that can be sampled grows exponentially. There is 
currently no consensus on the strategy of sampling this high-dimensional space. de Almeida Martins and 
Topgaard (2018)®° use a pseudorandom sampling strategy in their NMR experiments, but this is less trivial to 
set up on a clinical MRI scanner, where TE, TR, and TI are more commonly sampled on a grid2832,33,35 


Regardless of the sampling strategy, scan times for multidimensional experiments quickly become prohibitive. 
For efficient D-T2 NMR in homogenous media, Ahola et al. (2015)!0! propose to encode different diffusion 
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weightings along one spatial dimension by simultaneous application of a gradient pulse and the refocusing pulse, 
i.e. the strength of diffusion weighting varies linearly as a function of position. By combining this with a CPMG 
sequence, the 2D experiment can be performed in a single scan. However, in inhomogeneous samples such as 
biological tissue other strategies have to be used. Hutter et al. (2018)28 propose to acquire multiple readouts after 
a single excitation to perform efficient diffusion-T2* relaxometry, and efficient schemes to sample at multiple 
TT's (see also Chapter 9). 


In addition to scan time, hardware limitations further challenge diffusion-relaxation correlation experiments in 
humans. In diffusion-T2 experiments, the minimally achievable TE for a given diffusion weighting is generally 
large: on a typical clinical system the minimal TE is around 110 ms for b = 5000 s/mm? (Figure 7.5a). This does 
not leave much room to disentangle compartments with shorter T2. Using the latest ultra-strong gradient 
hardware102,103 the accessible b-TE space becomes much larger. Strong gradients can be leveraged to compute 
T2 of different compartments in human white matter>>!94, Combining strong gradients with a shorter image 
readout, diffusion of short T2 components such as myelin water might become measurable105,106 (Figure 7.5b). 


Figure 7.5. Increase of the accessible parameter space with ultra-strong gradients and shorter readouts. 


90° 180° 


300 mT/m 80mT/m 40mT/m 


40 60 80 100 120 140 
TE [ms] 


a). b-value vs minimum achievable echo time (TE) for different maximum gradient amplitudes of 40 mT/m, 80 mT/m and 
300 mT/m (represented by different colours), and different readouts (represented by different line styles: solid lines EPI readoud 
and dashed lines spiral readout. 


Figure created by Dr. Lars Mueller and partially adapted from195 with permission from the authors. 
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TE [ms] 


b). With ultra-strong gradients and spiral readout it becomes possible to achieve short TE and high b-values. 


106 


Figure adapted from “° with permission from the authors. 


7.5. How to estimate diffusion and relaxation properties? 


A heterogeneous material is often described as a sum of microscopic domains that do not exhibit exchange with 
one another (see Chapter 1). In this framework, each micro-environment is characterised by its diffusion 
properties (in the absence of restriction/intra-compartmental diffusional kurtosis or at sufficiently low diffusion 
weightings this can be represented by diffusion tensor D, see Chapter 2) and relaxation properties (represented 
by T1, T2, and T2*, ignoring the frequency offset). The distribution P(D,T1,T2,T2*) then gives the probability of 
finding a given D,T1,T2 and T2*, and as such characterises the composition of micro-environments within a 
voxel. The signal can in this case be written as: 


S(B, TR, TE, TD) = S(0,0,0,0 Ha Í [ Í 
Sym (3) 


K (B, TR, TE, TD, D, T1, T2, T2*) P(D, T1, T2, T2*) [6] 
dT1 dT2 dT2* dD. 


Here, S(0,0,0,0) is the signal when no diffusion- or relaxation-weighting is applied, and Sym*(3) denotes the 
manifold of 3 x 3 symmetric positive-definite matrices. 


The functional form of the kernel K depends on the pulse sequence and the sampling scheme (i.e. which 
parameters are varied). For the spin-echo experiments without inversion recovery as described in Section 7.4.1, 
one can derive 
TR TE TD 
x) = — R: — 
K(B, TR, TD, TE, D,T1,T2,T2*) = exp(— B: pJ: exp[- TĪ 77) exe(- T> T> )ex0(- +=) . [7] 
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In the equations above, we can recognise the multi-dimensional Laplace transform formulation which was first 
used in the context of NMR in the early 90’s*!07, and only became widely used a decade later11,12,108 This 
formulation was later adapted to accommodate the distribution of diffusion tensors rather than 
diffusivities109-113_ and extended to correlate diffusion tensors with relaxation®. 


7.5.1. The inverse problem 


Eq. 6 formulates the forward problem, i.e. given a probability distribution P the resulting signal can be 
computed. Instead, we are interested in solving the inverse problem: given a set of signals S,, (which can be 
stored as a as an N x 1 vector s with n = 1,...,N) obtained by varying the acquisition parameters, i.e. 
(B,TR,TE,TD),, find the distribution P. This turns out to be a very hard problem to solve for the formulation 
above: it requires taking the inverse Laplace transform. The inverse Laplace transform is well-known to be ill- 
posed and ill-conditioned™: there is no unique solution! !4, and small errors in the initial data can result in 
much larger errors in the estimates. In contrast, the Fourier transform is well-defined and more easily 
implemented. 


A rich literature exists on solving the inverse Laplace transform, and the approaches can roughly be subdivided 
according to the priors they impose!!>>! 16. on the one hand, nonparametric inversion approaches typically do 
not impose any functional form on P other than that it should be positive and between certain bounds. This can 
be further extended with priors that promote smoothness and/or sparsity. On the other hand, model-based 
approaches impose functional forms of P following from a theory about the material, e.g. the number of 
compartments. Chapters 11 and 12 review these approaches into detail, here they will be briefly introduced and 


their differences discussed. In addition, we discuss methods that aim to estimate the moments of P. 


7.5.1.1. Nonparametric inversion 


The nonparametric approach discretises the solution P by pre-specifying sets Dm, T1imy,T2m3, T2in, With my = 1, 


.. My; m = 1,...,M2; m3 = 1,...,M3; m4 = 1,...,M4. The values in these sets are chosen such that they are within 
plausible bounds for the solution. Then, it aims to estimate the probability of finding pre-specified values of 
Dmp Tlm, T2m, T2iny Le. the discrete fractions fin mymymg Corresponding to each of the combinations of 


Dos Plig T2 125. 


To illustrate this, let us first consider a 1D experiment, i.e. where only TE is varied to estimate T2, and B = 0, TD 
= 0, TR > T1. In human brain, we would expect the T2 to be somewhere between 10 and 150 ms, so we could 
for example choose [T21,T22,..., T215] = [10,20,...,150] ms. In many applications, however, these values are not 
chosen to be linearly spaced but instead to be equally spaced in 1n(R2)! 17 to cover a wide range of decades, 
without the need to make the number of values M very large. When the fractions fm (sometimes also referred to 
as the ‘spectral amplitudes’ or ‘weights’ corresponding to each of the T2,,) are stored into an M x 1 column 
vector f with m = 1,...,M, Eq. 6 can be written as: 

s=Kf+e. [7] 
Here, K is an N x M matrix with pre-computed values of the kernel, i.e. Kym = exp(- 7} K is often called the 
kernel matrix or dictionary. e is an N x 1 column vector representing the experimental error in each 
measurement (e.g. due to noise). Note that if S is not normalised by S(0), the fractions fm do not sum up to unity 
but to S(0). 


In a multidimensional experiment, the total number of a priori defined solutions, i.e. the fractions finjmym3mq 
corresponding to each of the combinations of Dn, Tlm, T2m, T2iny is Mj M2M3M4. For convenient notation, 


this could be vectorised! into an M x 1 column vector f with elements fm and m = 1,...,.M where M = 
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MıM2M3M4. Here, Knm = exp(—B,: D,)(! — exp(- = )xp[- 772 )exp(- ==) and the system of equations 


can thus be written as Eq. 7. Noticeably, the kernel K quickly becomes enormous. Therefore, nonparametric 
inversion has mostly been applied in 1D or 2D experiments. Some solutions to deal with large sizes of K have 
been proposed that either compress K in the measurement direction” or in the ‘signal predictions direction. 


Regarding compressing the data, a commonly used solution proposed by Venkataramanan and Song”! first 


consider the sub-kernels K; of each dimension in case they are separable (e.g. K> nm = 1 — exp(- T) and 

TE, ere 
K3 nm = exp(- T] for T1 and T2, with dimensions Nz x Mo and Ns x Ms respectively, where Nz and N3 
denote the number of unique values of TR and TE in the experiment). They then reduce the dimensions by 
singular value decomposition (SVD). Specifically, SVD is a generalisation of the eigendecomposition and 
factorises a real matrix as K = UXVT, where X is an N x M matrix with the singular values (non-negative real 
numbers) on the diagonal in descending order, and U and V are unitary matrices (i.e. UTU = UUT = I with I the 
identity matrix) and have size N x N and M x M respectively, with the columns containing the left- and right- 


singular vectors. By truncating X to 2 with size S x S keeping only a subset of the largest singular values, the 
reduced kernel of size S x M becomes K = EV", where V is also truncated to have size M x S. Doing this for 
each sub-kernel, the problem can finally be written as Eq. 7 with K = K> @ K; € aad 
the tensor product. s contains a compressed version of the data, with size S283 x 1. This is a vectorised version of 
S= 0! SU> ER? $3, where Ü; is truncated to have size N; x S; and S € RN2*N3 
following each combination of parameters TR and TE. 


where ® denotes 


contains the signal 


While the approach above requires the acquisition of each combination of parameters, a lot of data is 
immediately thrown away by the SVD compression! 18119, Acquiring all combinations is time-consuming, and 
they propose instead to acquire only a few combinations and complete the matrix ,Š by adopting compressive 
sensing techniques. Mitchell and Fordham (2011)!2° highlight difficulties with the sub-kernel SVD approach 
when the kernel K is not separable, and instead propose to compress the data by window averaging: at long TE, 
for example, all short T2 components have decayed away and the number of data points can be reduced by 
progressively averaging more data points as TE increases. 


Regarding reducing the possible number of solutions, de Almeida Martins and Topgaard (2016)° propose to 
sparsely sample the solution set and perturb the solutions iteratively, so that the size of the kernel in each 
iteration remains limited. For the diffusion tensors, which theoretically add 6 dimensions to the kernel, they 
additionally assume axial symmetry so that it can be described by its perpendicular and parallel diffusivity, and 
two angles defining the first eigenvector (D | ,Dj|,n(0,¢)). Yang et al (2018)!2! implement randomized SVD in 
which only a single solution at a time is computed to update the SVD process. In addition, they use more coarse 
dictionaries and obtain intermediate solutions by polynomial fitting. 


Once the problem is written in the form of Eq. 7, one could in principle get an estimate f by numerically 
computing the pseuso-inverse K*, for example from SVD as Kt = VEtU!. f can then be estimated as 
f=Kts. [8] 


1 Vectorisation allows a tensor T of order L and dimensions J4, J2, ..., JL with elements aj, ;,.._;, to be interpreted as a 


T z š 
vector such that t= [ai1.1.21...1 TRA QT 1.1 AJ12...1> es aJ|J2...J L| = [ų,m,...az] with J = J4, J2, ..., JL, i.e. 


lexicographically ordering the elements. 
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This formulation solves the problem in the least-squares sense: it is equivalent to minimising the sum of squared 
F 2 I : š I I 
residuals f = arg min z| |s -Kf | E It can be shown that the least-squares estimator is the optimal estimator (in 


the sense that it has the minimum variance of all estimators) if the experimental errors are uncorrelated and 
have zero mean and constant variance. For the problem at hand, however, this strategy will not give a very useful 
answer. The matrix K is ill-conditioned, which means that the inversion process will amplify the errors in S to 
large errors in f. In addition, there exist infinitely many solutions because not each of the columns of K are 
linearly independent. This can only be ameliorated by making some prior assumptions about f (i.e. 
regularisation), without necessarily imposing a functional form as in model-based approaches. 


Lawson and Hanson (1995)122 imposed a positivity constraint which increases the stability of the solution and 


coin their algorithm non-negative LS (NNLS). Provencher (1982)123 in addition proposed to add a 
regularisation term 


Î = agmini lls- KÈ + asf). 9] 


with g a penalty function and ç a regularisation parameter that can be tuned. Several works proposed to 
promote parsimony by employing so-called Tikhonov regularisation: a smooth solution with a minimum number 
of peaks is favoured. In practise, this means that a term is added to the objective function that is a function of the 


L2-norm off, i.e. g(f) = || f| a with T some suitably chosen Thikonov matrix. Daducci et al. (2015)124 in 


115,123 


addition add the L1-norm of f, i.e. || f |Ẹ multiplied by a second regularisation parameter $. The L1-norm further 
promotes sparsity of f, i.e. f only having a few non-zero elements, compared to the L2-norm!2°, and L1 or L2 
regularisation can be adopted depending on the application*®»!2®, Although regularisation makes the 
inversion less ill-posed, it might cause bias!2”. Setting the regularisation parameters is generally not trivial!28; it 
can either be set empirically or based on simulations, or derived from techniques such as generalised cross- 
validation??? or L-curve!39, 


Other ways to stabilise the solution have been proposed. Benjamini and Basser (2016)* suggest marginal 
distributions constrained optimisation (MADCO) in which they first estimate the 1D distributions, and then use 
these to constrain the estimation of the 2D spectrum. They showed that this approach could greatly reduce the 
amount of data needed. It remains to be evaluated how instabilities and degeneracies in the single-contrast 
solutions affect this result (see Section 7.3.3). Sun and Dunn (2005) 131 similarly constrained 3D inversion by 
first solving a 2D inversion (e.g. T1-T2), and then use the resolved spectrum to constrain the remaining 
inversion (e.g. D-T2). Kim et al. (2017)®! use the observation that in an imaging experiment, the spectra in 
neighbouring voxels should be correlated. They adopt spatial regularisation and jointly estimate the parameters 
for multiple voxels simultaneously rather than for each voxel individually. Alternatives to least-squares 
minimisation have also be proposed: e.g. maximum entropy inversion132,133_ de Almeida Martins et al.® note 
that including several B-tensor shapes in the acquisition (as opposed to the conventional Stejskal-Tanner 
experiment) stabilises the solution. 


7.5.1.2. Model-based estimation 


A theory about which degrees of freedom have the most relevant contribution to the signal allows one to develop 
a model: a concrete realisation of the theory that describes how the signal is related to the relevant degrees of 
freedom!54, In contrast, representations can be regarded a mathematical means to merely compress or describe 
the signal. Following this definition, from the inversion procedure discussed in the previous section we obtain a 
representation of the signal; a distribution of Gaussian signal-contributions. Even if a functional form of this 
distribution were to be imposed, this would not necessarily be classified as a model according to the definition 
above if a clear theory is lacking; it might simply be out of convenience. 
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Here we take a more pragmatic approach and discuss several models and representations that impose a 
functional form of the distribution P. The continuous distribution fit has shown to be very sensitive to the 
SNR!0135, and introducing further priors can significantly stabilise the fit. For example, if the peaks are narrow, 
imposing delta priors describes the data as well as the continuous distribution, but with fewer parameters. Eqs. 6 
and 7 have been at the heart of several techniques in human brain white matter (in the absence of intra- 
compartmental diffusional kurtosis and exchange); Table 7.1 gives an overview. 


Several priors have been proposed for T2 fitting. Does and Gore (2000)®° consider an SNR per pixel of 
approximately 30 at the first echo point impractical for nonparametric inversion, and impose delta priors to fit a 
one-, two-, and three-component slow-exchange model to the data, where an F-test was used to evaluate which 
model best described each dataset. Stanisz and Henkelman (1998)>4 instead impose a Gaussian prior on the 
distribution of T2 on the logarithmic scale, because delta priors did not give a statistically correct result. 


For diffusion data, a wide range of techniques have been proposed that formulate the signal as resulting from a 
sum of microenvironments each represented by a tensor136, and we will discuss a few examples here. Tuch et al. 
(2002)!37 represent crossing white matter fascicles in the brain by different tensors and fix the diffusivities of the 
tensors, while they estimate the direction of each fascicle. Behrens et al. (2003)!38 estimate two tensors per voxel; 
one ‘stick for the axonal compartment with zero perpendicular diffusivity and one ‘ball’ with similar 
perpendicular and parallel diffusivity. Zhang et al. (2012)!3? estimate three diffusion tensors (a ‘stick for the 
intra-axonal compartment, ‘zeppelin’ with unequal perpendicular and parallel diffusivity for the extra-axonal 
compartment, and one ‘ball’ for cerebrospinal fluid), and assume that the first eigenvectors are distributed 
according to a Watson distribution. They fix the parallel diffusivities and employ a turtuosity assumption for the 
extra-axonal compartment!40, Jelescu et al. (2016)8! investigate the stability of model parameter estimation for 
the model in Zhang et al. (2012)!39 without the ‘ball’ compartment, and with free estimation of all the 
diffusivities. They find a model degeneracy when fitting this model up to moderate b-values; two biophysically 
plausible solutions fit the data equally well. 


Instead of imposing delta priors on diffusivities as in many of the examples described so far, other distributions 
have been proposed. Lasič et al. (2014)!4! first compute the so-called powder-average, which is the averaged 
signal over all directions, to factor out the dependence of the signal on orientation dispersion. They then propose 
to model P(D) with a Gamma distribution 142143, Scherrer et al. (2016) !13 represent each fascicle by a 
distribution of tensors, which they parameterise by a matrix-variate Gamma distribution. Interestingly, the 
concentration parameter in the Gamma distribution can capture the overall compartment’s heterogeneity arising 
from different sources, e.g. orientation dispersion or heterogeneity in underlying diffusivities. 


For combined diffusion-relaxation correlation data, early works first try to find evidence for multiple 
components from a single-contrast, and subsequently obtain estimates for each component in the second 
dimension. Van Dusschoten et al. (1995)10 estimate the number m of ensembles, and the signal fraction and 
corresponding T2 from a CPMG T2 measurement. They then perform n > m PGSE experiments with n different 
TEs, and estimate n Dapp values. For each PGSE experiment, the relative amplitude of the m fractions can be 
calculated given the TE of that experiment and the T2 corresponding to each fraction. Obtaining estimates of 
the diffusivities for each ensemble then comes down to solving n linear equations with e.g. a SVD routine. 
Beaulieu et al. (1998)°> estimates D per TE from a CPMG prepped PGSE sequence, approximating the signal as 
a monoexponential decay in b. 


For the correlation of diffusion tensors and relaxation values, De Santis et al. (2016)32 first estimate the number 
and directions of white matter fascicles in each voxel using spherical deconvolution with an a priori defined 
kernel!5151_ They then fix the directions and estimate the T1 and diffusivities of the tensor associated with 
each fascicle, where they further reduce the number of parameters by using the tortuosity assumption. Veraart et 
al. (2018)? reduce the dimensionality of the parameter space by first computing the powder average of the 
diffusion data per TE, and then jointly fitting a ‘stick-zeppelin’ with two different T2’s to the diffusion and 
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relaxation data. As in previous work®! 


potential to resolve the degeneracy. 


they do not fix diffusivities, and show that adding the T2 dimension has 
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7.5.1.3. Estimation of moments 


The free inversion and model-based estimation approaches typically aim at characterising the distribution P 
either by imposing some properties (e.g. smoothness) or a functional form. Instead of characterising the 
distribution itself, several studies have proposed to estimate its moments (or cumulants, which provide an 
alternative to the moments) as useful descriptors of its shape!>2. The cumulants can be directly estimated from 
the signal without the need to first estimate P. 


In diffusion MRI this strategy is commonly adopted!53-!54, The cumulants are related to the Taylor expansion in 
b of the logarithm of the signal near b = 0, i.e. in 1D: 


InS(b) = InS(0) — bD — PPD K +- [10] 


Here D is the mean of the distribution, and K is the excess kurtosis, related to the variance of the distribution155 


p= f prpup [11] 
0 


var(D) = f “(D- Dy P(D)dD 2] 


Diffusion tensor imaging (DTI!4®) considers only the first two terms in Eq. 10, and diffusion kurtosis imaging 
(DKI!>>) involves the first three terms. It should be noted that the Taylor expansion does not converge for high b 
in the case of multiple Gaussian compartments (Eq. 7)!>3; the accuracy of the estimates will decrease when the 
b-values are beyond the radius of convergence!®. 


Estimation of moments or cumulants has been performed in relaxometry!°” and multidimensional diffusion 
MRI!!2 separately, and very recently also for D-T2 relaxometry data”?. 


7.5.2. Uncertainty 


Most works report a single solution to the inverse problem; this point estimate aims to reflect the “best guess” of 
the parameter value. Point estimates, however, do not provide any information on the uncertainty or 
interdependency of parameters. For nonparametric inversion, Prange and Song (2009)!°8 aim to capture the 
uncertainty and estimate an ensemble of solutions using Monte Carlo sampling. Specifically, based on Bayes’ 
rule, they sample from the posterior distribution P(f|s) < P(s|f)P(f), with P(s|f) the likelihood function and P(f) 
the prior. The likelihood function contains a formulation for the forward problem, e.g. s = Kf as in Eq. 7, anda 
noise model. Markov Chain Monte Carlo (MCMC) algorithms are commonly used to sample from the posterior 
distribution. Prange and Song (2009)!°8 propose a non-negative multi-normal distribution for the likelihood 
and a constant prior giving an equal probability to any value of f. Upon sampling the posterior distribution of a 
T2 experiment with an enhanced Gibbs sampler!>9, they find an increasingly upward bias in spectral values with 
decreasing T2, and correlations between spectral positions. Several studies?>-62,69.85 use bootstrap procedures 
which randomly sub-select data points and compute a spectrum for each subsampling. From this ensemble of 
solutions, confidence intervals can be estimated and Laplace inversion can be stabilised against spurious peaks in 
the correlation maps. 


For model-based estimation, the principles above have similarly been used to quantify uncertainty and 
correlations, e.g. using MCMC to sample from the posterior distribution of several microstructural 
models!38:160-162 Harms and Roebroeck (2018)! provide an efficient GPU-based implementation of MCMC 
that is compatible with a wide range of models. Alternatively, several works164-166 use bootstrapping to 
characterise uncertainty in diffusion MRI estimation. Such estimates of uncertainty have been integrated in 
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tractography for the virtual reconstruction of fibre pathways: probabilistic tractography - as opposed to 
deterministic tractography - aims to take into account the uncertainty of the fiber direction estimate in each 
voxel. 


7.5.3. Laplace vs Fourier 


In addition to the Laplace transform in Eq. 6, the Fourier transform is used for multi-dimensional experiments 
involving diffusion as well; an excellent overview is given in Callaghan (2011)!. In a 1D diffusion NMR 
experiment with narrow pulses, the Fourier transformation with respect to q (with q = yóG where y is the 
gyromagnetic ratio, ô is the pulse duration, and G is the gradient strength) gives the spectrum of 
displacements!®’>168_ For Gaussian diffusion with diffusion coefficient D, this spectrum is a Gaussian function 
with width D. In the case of multiple microenvironments, each environment i would exhibit a Gaussian 
spectrum with width Dj, and the overall spectrum would be a superposition. This spectrum would be very hard 
to disentangle, and as such the spectrum of diffusion coefficients P(D) (or P(D) in the case of the distribution of 
tensors in Eq. 6) might be more intuitive to work with. 


Several studies report combined Laplace- and Fourier transformation in diffusion-relaxation correlation. Britton 
et al. (2004)169 correlate T2 with molecular displacement in water diffusing through an alginate bead pack. They 
performed Fourier transformation to obtain a propagator for each TE, and inverse Laplace transformation for 
each displacement in the propagator. Washburn and Callaghan (2007)!79 perform diffusion T2-T2 exchange 
experiments and combine a Fourier dimension with two inverse Laplace dimensions. As such, they can better 
determine whether changes in T2 comes from intra- or inter-pore transport, as they have access to the 
displacement of spins. 


7.6. How to validate diffusion-relaxation correlation? 


Eq. 6 in the previous section was presented as an overarching framework of methods that describe materials as a 
sum of microscopic domains that are each already fully coarse grained, which means that the diffusion time is 
long enough for the water molecules to have already homogenised intra-compartmental microstructure, and are 
not in exchange. If these assumptions are invalidated, estimates from inversion- or model-based analysis could 
become biased and/or meaningless. In this Section, we briefly discuss strategies that have been used to validate 
and evaluate results from diffusion-relaxation measurements (for a more general overview see e.g. 171-172 and 
references therein), and hallmarks that could point to invalid assumptions. Specifically, we will discuss 1) 
computer simulations, 2) phantoms, 3) histology, 4) adjusting the measurement regime, and 5) adding MRI 
dimensions. Chapters 13 and 14 consider validation into more detail. 


7.6.1. Computer simulations 


With the development of improved computational hardware, computer simulations have become an important 
tool to predict and validate MR measurements. Assuming a ground truth, forward simulation of Eq. 6 is 
commonly used to evaluate precision and accuracy of the inversion®?»!73, For more complex scenarios, such as 
the presence of restrictions, various approaches have been developed. Some rely on the evaluation of Fick’s 
second law of diffusion; given a geometry, this partial differential equation can be solved analytically for a 
limited range of initial- and boundary conditions, or numerically for more general configurations. Brownstein 
and Tarr (1979)!74 and Callaghan (1995)!79 derived expressions for relaxation- and PGSE experiments 
respectively, in planar, cylindrical, and spherical pores under conditions of surface relaxation based on the 
eigenmode expansion of the propagator. This was also extended to general waveforms by breaking the waveform 
into intervals and writing a propagator for each interval!76-!79, Using this framework, Callaghan et al. 
(2003)!89 simulated T2-D signals from restricted geometries (plane and sphere) with varying size and wall 
relaxivities, and study the 2D Laplace inversion spectra’. The spectra showed, in addition to the primary mode, a 
wide spread of diffusion and relaxation values even for these simple pores (Figure 7.6a)!8!. These extra features 
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might prove of value to separate surface- from bulk fluid in the pore. Novikov et al. (1998)!82 developed a 1D 
numerical model of restricted diffusion and relaxation in multi-compartment planar geometries for PG CPMG, 


and this was later extended to cylindrical geometries! ®3, 


Other approaches perform Monte Carlo simulations of the movements of a large number of hypothetical 
molecules. These approaches are becoming increasingly popular because of their greater flexibility to model 
arbitrary configurations and multiple processes simultaneously (e.g. diffusion, exchange, susceptibility, and 
relaxation). For example, Lin et al. (2018)”° simulate myelinated axons with different T2 in each compartment, 
and diffusing molecules can exchange between compartments with a probability according to a pre-defined 
axonal membrane permeability (Figure 7.6b). Ruh et al. (2018)°° use Monte Carlo simulations to study diffusion 
in microenvironments with heterogeneous susceptibility. 


Figure 7.6. Computer simulations in diffusion-relaxation correlation. 
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a). Simulation of a simple planar restricted geometry with wall relaxation gives complex features in a D-T2 plot. The diagonal 
arrow indicates the position of the primary relaxation—diffusion mode obtained from the low-q data. 


Figure adapted from!®° with permission from Elsevier, Copyright 2003. 
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b). Monte Carlo simulation geometry of myelinated axons in a periodic array of hexagonal columns, each axon is surrounded by a 
permeable myelin sheath. 


Figure adapted from” with permission from John Wiley and Sons, Copyright 2018. 


7.6.2. Phantoms 


Phantoms to validate diffusion-relaxation correlation have a well-characterised composition and are usually 
man-made and non-biological, although food phantoms have also been used (Chapter 13). Realistic phantoms 
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typically contain an NMR visible liquid with similar diffusion coefficient and relaxation times as can be found in 
the material of interest, and an NMR invisible material that reflects realistic architectural configurations, e.g. 
intra- and extra-axonal spaces and/or crossing fibres for white matter. Regarding NMR-visible liquids, several 
examples can be found in literature to validate multi-dimensional MR experiments, e.g. to evaluate strategies for 
faster acquisition and analysis against conventional methods. Benjamini and Basser (2016, 2018)1:126 used 
doped water and polyvinylpyrrolidone (PVP) to create T1-D phantoms with distinct peaks, which were 
characterised individually in advance. Ahola et al. (2015)!9 use porous silica gel powder immersed in water to 
create D-T2 spectra from water in the pores and between the particles. Dairy cream has been used as a phantom 
that has both bi-exponential diffusion as well as T2!6184, 


Regarding the architectural configuration, careful thought should be given to the dimensions, surface relaxation, 
and susceptibility effects of the material in which the liquid will reside. de Almeida Martins and Topgaard 
(2018)85 combine a liquid crystal exhibiting anisotropy!®> with a yeast suspension composed of microscopic 
spherical cells with intra- and extra-cellular isotropic diffusion of different magnitudes. They use this to evaluate 
6-dimensional T1-T2-D inversion (Figure 7.7). 


7.6.3. Histology 


Histology can be used to study the microscopic anatomy of biological tissues and thus provides an independent 
window into the composition. Tissues are typically fixed to preserve their structure and then embedded in a 
medium. Sections of the tissue can then be stained to give contrast to the tissue and highlight specific structures 
or chemical components. Benjamini and Basser (2017)? used immunohistochemistry to mark astrocytes, 
microglia, myelin oligodendrocyte glycoprotein, and neurofilaments with antibodies, in sections of a rat spinal 
cord. They then correlated properties of T2-D derived components to the relevant immunohistochemistry 
features. Although this cannot be used as an absolute “ground truth” for MRI-derived measures because of 
differences in the source of the signal and changes after fixation and embedding, it can provide useful 
independent information (Chapter 14). 


7.6.4. Adjusting the measurement regime 


Modelling all the compartments in a material typically involves the estimation of a relatively large number of 
parameters. If prior knowledge exists on differences in diffusion- and/or relaxation properties between 
components, the experiment can be set up to ‘filter out’ certain compartments while keeping others. This might 
facilitate the fitting and rule out errors from the inversion procedure; the results from the full- and filtered 
experiment should at least be compatible. Sorland et al. (2004)13 used this principle in a T2-D experiment of 
crude oil and water in rock cores, where they apply sufficiently high diffusion gradients to filter out the signal 
from water and study the T2 decay of oil. McKinnon and Jensen (2019)3! and Kleban et al. (2019)76 used high 
diffusion-weightings to sufficiently decay the signal from extra-axonal water and study the T2 and T2* from 
intra-axonal water, respectively. Tax et al.186 used high b-values and isotropic diffusion weighting to filter out 
compartments in which water can move at least in one direction, thereby targeting small spherical 
compartments such as cell bodies. Andrews et al. (2006)°° targeted myelin water by using an inversion pulse 
with suitable TI to nullify signal from the intra- and extra-axonal spaces in excised frog sciatic nerve. They then 
measure diffusion and T2 from the remaining compartment, and compare this to measurements that are 
sensitive to all compartments simultaneously. 


7.6.5. Adding MRI dimensions 


Apart from diffusion and relaxation, other phenomena play a role at the time scale of typical measurements. For 
example, exchange between compartments can affect the observed fractions!8’, and mesoscopic field gradients 
can bias the estimated relaxation rates. The influence of these effects might be estimated with additional MRI 
measurements, effectively adding dimensions beyond diffusion and relaxation. To measure exchange, the same 
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Figure 7.7. Phantoms for diffusion-relaxation correlation. This phantom is composed of a yeast suspension microstructure (green 
frame, set of spherical shells representing the cell membranes that delimit the intra- and extracellular domains), and a liquid crystal 
microdomain (red frame, coherent set of cylindrical water channels in a continuous matrix of detergent and hydrocarbon. The blue 


27 


frame demarks the entire sample volume. Six dimensional P(R1, R2, Diso» D||/D 1, 9, b) distribution calculated through unconstrained 


inversion of data acquired with the phantom. 
Copyright 2018 


Figure adapted from” under the Creative Commons license (http://creativecommons.org/licenses/by/4.0/) 


property (e.g. relaxation or diffusion) has to be encoded at two different times188,189 If exchange takes place 
during the mixing time, the spectra will reveal peaks located off the diagonal. Chapter 5 considers exchange 
experiments into more detail. 


Regarding local internal magnetic field gradients, Carr and Purcell (1954)*? showed that the sensitivity of the 
measurement to the additional dephasing is dependent on the time between echoes in a CPMG sequence; 
shorter durations result in less echo attenuation. By incorporating a dimension in which this timing is varied, 
the local field can potentially be estimated. Seland et al. (2004)18 perform D-T2 experiments with an extra 
dimension for internal gradients in liquid-saturated packings of glass beads and perform Laplace inversion to 


disentangle the different phenomena. 


7.7. Summary and Discussion 


This Chapter has focused on combining diffusion and relaxation MR to study the heterogeneity of biological 
tissue. It should be apparent that relaxometry provides rich information that can complement the information 


obtained from studying the diffusion process. To date this has been mostly applied in excised tissue and animal 
imaging, but more recently diffusion-relaxation correlation has also found its way into in-vivo human imaging. 


Improvements in MRI hardware and imaging strategies, such as strong gradients, multi-band imaging, and 


efficient image readouts, have greatly facilitated this translation, but challenges remain that currently hamper its 


widespread adoption. In the next Sections, we will discuss some of the challenges and limitations, and provide an 


outlook for the use of diffusion-relaxation correlation in clinical applications. 
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7.7.1. Challenges and limitations 


A first limitation is noise. The signal should be measured until it has almost completely decayed, ideally by at 
least three orders of magnitude before it reaches the noise floor128 This would correspond to an SNR of at least 
1000, which is much higher than what is commonly achieved in in-vivo human imaging. This might to some 
extent be addressed by averaging multiple acquisitions, but this is in turn challenged by the limited acquisition 
time. Denoising techniques have been developed for MR images and for diffusion MRI specifically (e.g. 199191), 
but care should be taken that no relevant information is removed from the images. In addition to the presence of 
noise, the distribution of the noise can further challenge accurate estimation. The noise in MRI magnitude 
images are not Gaussian distributed; they can be Rician distributed or non-central Chi distributed, with 
stationary or non-stationary noise, depending on the combination of the coils and parallel imaging techniques 
used!92, The commonly used least-squares formulation can give biased results in this case. Nevertheless, 
promising results have been shown in in-vivo human data where the problem was sufficiently constrained, e.g. in 
model-based approaches27;:106 in 2D inversion approaches”®, and in higher dimensional inversion approaches 
with multidimensional diffusion33 


Establishing reasonable assumptions is another major challenge. Model-based approaches commonly rely on a 
theory (i.e. a ‘picture’ of what the reality might look like). Based on this, assumptions on the relevant degrees of 
freedom might be derived, for example diffusion in white matter is commonly modelled with two 
compartments. However, if such a theory is not available (e.g. in disease), an interplay between modelling and 
nonparametric inversion could provide a useful strategy to find suitable constraints. As inversion strategies have 
been so commonly used in NMR of porous media over the years, common artefacts are well-recognised and 
hallmarks of invalid assumptions can sometimes be derived from the spectra. For example, ‘pearling’ in which 
local islands of peaks arise, is well-known to be an artefact of the transformation!. An example of a hallmark is 
when both T1 and T2 are part of the multi-dimensional experiment, exchange can dislocate peaks and add 
additional ones that may have negative amplitudes incompatible with the positivity constraint!3, Restricted 
diffusion and wall-relaxation may also exhibit specific hallmarks!®°, Once suitable constraints have been 
established, a more accurate quantification could be obtained by fitting the compartmental model rather than 
reading the values from the spectra124. 


Perhaps the most challenging aspect of translating diffusion-relaxation correlation to in-vivo human imaging is 
the acquisition and the long measurement times associated with it. MR fingerprinting uses very fast acquisitions 
and relies on a large amount of lower quality images to derive the fingerprint?!. Other ways of fast acquisition 
are discussed in Chapter 9. Alternatively, combining compressive sensing techniques with sparse acquisitions 
can reduce the number of necessary data points. Finally, some practical challenges may complicate acquisition; 
commercially available diffusion sequences do not always allow to change TE and diffusion time independently, 
which unintentionally mixes these experimental dimension2’. Also setting up acquisitions with many different 
TEs and TRs can be cumbersome and require careful bookkeeping. 


7.7.2. Outlook 


Diffusion-relaxometry has the potential to provide a more complete picture of tissue microstructure than their 
1D counterparts. In clinical applications, the increased sensitivity and specificity to changes in tissue 
microstructure can result in improved disease monitoring and diagnosis. This links in with the trend of 
personalised medicine, in which the treatment is tailored to the specific characteristics of each individual. 
Preliminary applications in different organs, such as brain, prostate, and placenta, have shown that the potential 
usefulness of this technique is widespread. Moving forward, diffusion-relaxation correlation can further be 
extended in terms of dimensions; it is well-known that T2* in white matter is anisotropic and thus depends on 
the orientation of the tissue with respect to the main magnetic field, and recent studies also show T2 
anisotropy? t195, If different rotations of the tissue in the scanner are feasible, T2*- and T2 tensors can thus 
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potentially be estimated. Finally, data-driven and machine learning methods provide exciting new avenues to 


analyse the high-dimensional data in diffusion-relaxation correlation 


196, 
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